Differential gene expression

Material

Exercises

Find all markers for each cluster

Load the seu_int dataset you have created yesterday:

seu_int <- readRDS("seu_day2-4.rds")

And load the following packages (install them if they are missing):

library(Seurat)
library(edgeR) # BiocManager::install("edgeR")
library(limma)
library(dplyr)
library(scuttle)

The function FindAllMarkers performs a Wilcoxon plot to determine the genes differentially expressed between each cluster and the rest of the cells. Other types of tests than the Wilcoxon test are available. Check it out by running ?Seurat::FindAllMarkers.

Now run analysis:

de_genes <- Seurat::FindAllMarkers(seu_int,  min.pct = 0.25,
                                   only.pos = TRUE)
Time for coffee

This takes a while. Have a break.

Subset the table to only keep the significant genes, and you can save it as a csv file if you wish to explore it further. Then extract the top 3 markers per cluster:

de_genes <- subset(de_genes, de_genes$p_val_adj < 0.05)
write.csv(de_genes,
          "de_genes_FindAllMarkers.csv",
          row.names = F, quote = F)

top_specific_markers <- de_genes %>%
  group_by(cluster) %>%
  top_n(3, avg_log2FC)

And generate e.g. a dotplot:

dittoSeq::dittoDotPlot(seu_int,
                       vars = unique(top_specific_markers$gene), 
                       group.by = "integrated_snn_res.0.3")

Exercise

What are significant marker genes in cluster 0 and 8? Are the T cell genes in there?

You can re-load the vector with immune genes with:

tcell_genes <- c("IL7R", "LTB", "TRAC", "CD3D")
de_genes[de_genes$gene %in% tcell_genes,] |> knitr::kable()
p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
CD3D 0 2.0608957 0.769 0.227 0 0 CD3D
TRAC 0 1.6978332 0.619 0.204 0 0 TRAC
LTB 0 1.5374343 0.759 0.394 0 0 LTB
IL7R 0 1.5292274 0.437 0.114 0 0 IL7R
LTB1 0 1.1510176 0.673 0.465 0 7 LTB
TRAC1 0 1.8345938 0.759 0.274 0 8 TRAC
CD3D1 0 1.6980383 0.798 0.326 0 8 CD3D
LTB2 0 1.6381762 0.768 0.462 0 8 LTB
IL7R1 0 1.1453364 0.461 0.173 0 8 IL7R
LTB3 0 0.8247275 0.753 0.467 0 11 LTB

So, yes, the T-cell genes are highly significant markers for cluster 0 and 8.

Differential expression between groups of cells

The FindMarkers function allows to test for differential gene expression analysis specifically between 2 groups of cells, i.e. perform pairwise comparisons, eg between cells of cluster 0 vs cluster 2, or between cells annotated as T-cells and B-cells.

First we can set the default cell identity to the cell types defined by SingleR:

seu_int <- Seurat::SetIdent(seu_int, value = "SingleR_annot")

Run the differential gene expression analysis and subset the table to keep the significant genes:

deg_cd8_cd4 <- Seurat::FindMarkers(seu_int,
                                   ident.1 = "CD8+ T cells",
                                   ident.2 = "CD4+ T cells",
                                   group.by = seu_int$SingleR_annot,
                                   test.use = "wilcox")
deg_cd8_cd4 <- subset(deg_cd8_cd4, deg_cd8_cd4$p_val_adj<0.05)
Exercise

Are CD8A, CD8B and CD4 in there? What does the sign (i.e. positive or negative) mean in the log fold change values? Are they according to the CD8+ and CD4+ annotations? Check your answer by generating a violin plot of a top differentially expressed gene.

You can check out the results with:

View(deg_cd8_cd4)
p_val avg_log2FC pct.1 pct.2 p_val_adj
CD8A 0.0e+00 1.2956354 0.344 0.008 0.0000000
CTSW 0.0e+00 0.9686281 0.283 0.029 0.0000000
CCL5 0.0e+00 2.5884630 0.292 0.061 0.0000000
CD8B 0.0e+00 0.8536693 0.479 0.177 0.0000000
NKG7 0.0e+00 2.0692198 0.231 0.037 0.0000000
CST7 0.0e+00 0.9061154 0.147 0.012 0.0000000
GZMA 0.0e+00 1.2771487 0.173 0.025 0.0000000
TRGC2 0.0e+00 0.9457918 0.148 0.017 0.0000000
ID2 0.0e+00 0.6943385 0.579 0.347 0.0000000
MT-CO1 0.0e+00 0.3476088 0.988 0.980 0.0000000
KLRD1 0.0e+00 0.7547626 0.110 0.008 0.0000000
HCST 0.0e+00 0.6611536 0.694 0.489 0.0000000
GZMK 0.0e+00 0.8785594 0.134 0.022 0.0000000
TRGC1 0.0e+00 0.7533633 0.113 0.013 0.0000000
MT-ND4 0.0e+00 0.3442778 0.965 0.931 0.0000000
RP11-291B21.2 0.0e+00 0.5729470 0.226 0.076 0.0000000
MT-CO2 0.0e+00 0.2568810 0.994 0.994 0.0000000
FHIT 0.0e+00 -0.6810903 0.113 0.268 0.0000000
TRDC 0.0e+00 0.9507650 0.107 0.018 0.0000000
LYAR 0.0e+00 0.6692540 0.203 0.074 0.0000000
ACTB 0.0e+00 0.3240542 0.965 0.926 0.0000000
PRF1 0.0e+00 0.4148957 0.102 0.017 0.0000000
MT-CO3 0.0e+00 0.2992841 0.971 0.950 0.0000000
CD4 0.0e+00 -0.4000835 0.012 0.103 0.0000000
CCL4 0.0e+00 1.1567053 0.120 0.030 0.0000000
MT-ND2 0.0e+00 0.3127492 0.961 0.916 0.0000000
AC092580.4 0.0e+00 0.5206621 0.174 0.062 0.0000000
HLA-B 0.0e+00 0.3221429 0.974 0.932 0.0000001
RUNX3 0.0e+00 0.3428665 0.184 0.070 0.0000001
TPST2 0.0e+00 0.4343363 0.141 0.044 0.0000001
NR4A2 0.0e+00 0.5202848 0.350 0.199 0.0000003
MT-ND3 0.0e+00 0.2809915 0.958 0.928 0.0000005
GNLY 0.0e+00 1.6651734 0.135 0.045 0.0000006
MT-ATP6 0.0e+00 0.3132901 0.938 0.899 0.0000015
FAM173A 0.0e+00 0.4097859 0.197 0.086 0.0000022
CBLB 0.0e+00 0.2888085 0.131 0.044 0.0000091
BZW1 0.0e+00 0.4436429 0.360 0.221 0.0000109
JUN 0.0e+00 0.5248833 0.691 0.582 0.0000262
IL32 0.0e+00 0.4441875 0.742 0.624 0.0000293
LITAF 0.0e+00 0.4384073 0.417 0.272 0.0000299
IFRD1 0.0e+00 0.4287728 0.311 0.179 0.0000357
ACTG1 0.0e+00 0.3550782 0.750 0.619 0.0000501
LCP1 0.0e+00 0.3834983 0.318 0.196 0.0002808
HLA-DPB1 0.0e+00 0.3699177 0.215 0.111 0.0003358
KLRG1 0.0e+00 0.4372144 0.106 0.037 0.0003461
CYBA 0.0e+00 0.3626436 0.692 0.575 0.0005623
RPL36A 0.0e+00 -0.2636701 0.988 0.990 0.0006252
DUSP1 0.0e+00 0.4773763 0.700 0.576 0.0007094
DUSP2 0.0e+00 0.6941207 0.325 0.210 0.0008438
S100A4 1.0e-07 0.4158906 0.591 0.467 0.0009771
SRGN 1.0e-07 0.5274908 0.476 0.355 0.0010323
MT-ND5 1.0e-07 0.3723932 0.729 0.644 0.0012021
GZMM 1.0e-07 0.3179818 0.360 0.229 0.0013519
MT-ND1 1.0e-07 0.2968126 0.849 0.787 0.0015110
ZFP36 1.0e-07 0.4700353 0.399 0.274 0.0016737
TRAT1 1.0e-07 -0.5086304 0.138 0.238 0.0016743
MAL 1.0e-07 -0.4840650 0.152 0.252 0.0021382
A1BG 1.0e-07 0.3035121 0.149 0.068 0.0026934
ARPC5L 2.0e-07 0.2581570 0.176 0.086 0.0029573
HLA-A 2.0e-07 0.2922902 0.898 0.841 0.0032584
ICOS 2.0e-07 -0.3882220 0.036 0.105 0.0033535
S100B 2.0e-07 0.5873272 0.104 0.040 0.0046322
LINC00152 3.0e-07 0.3001131 0.104 0.040 0.0055603
STK17A 3.0e-07 0.3529199 0.430 0.306 0.0057236
C12orf75 3.0e-07 0.3524864 0.202 0.111 0.0064928
TSPAN32 4.0e-07 0.3176836 0.125 0.054 0.0068297
TSPYL2 4.0e-07 0.3144523 0.274 0.166 0.0071097
ARPC1B 4.0e-07 0.3041903 0.530 0.404 0.0078444
SH3BGRL3 5.0e-07 0.2899331 0.720 0.616 0.0089486
CXCR4 5.0e-07 0.3576632 0.739 0.634 0.0089842
RHOA 8.0e-07 0.3034162 0.427 0.297 0.0140834
CMC1 8.0e-07 0.4137942 0.103 0.041 0.0148046
SRSF7 1.0e-06 0.3534034 0.553 0.434 0.0178885
PSME2 1.2e-06 0.3612496 0.437 0.318 0.0227285
RAB27A 1.5e-06 0.2514754 0.109 0.045 0.0281147
MYO1F 1.5e-06 0.3325859 0.149 0.074 0.0287064
DNAJA1 1.8e-06 0.3733802 0.309 0.206 0.0338759
ADGRE5 1.9e-06 0.2649811 0.203 0.114 0.0346912
CORO1B 2.6e-06 -0.4747401 0.125 0.206 0.0477291

For an explanation of the log fold change have a look at ?Seurat::FindMarkers. At Value it says:

avg_logFC: log fold-chage of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group

To view CD8A, CD8B and CD4:

deg_cd8_cd4[c("CD4", "CD8A", "CD8B"),]
            p_val avg_log2FC pct.1 pct.2    p_val_adj
CD4  1.070126e-13 -0.4000835 0.012 0.103 1.998246e-09
CD8A 1.409306e-77  1.2956354 0.344 0.008 2.631597e-73
CD8B 7.113148e-36  0.8536693 0.479 0.177 1.328238e-31

Indeed, because we compared ident.1 = “CD8+ T cells” to ident.2 = “CD4+ T cells”, a negative log2FC for the CD4 gene indicates a lower expression in CD8+ T-cells than in CD4+ T-cells, while a positive log2FC for the CD8A and CD8B genes indicates a higher expression in CD8+ T-cells.

Plotting the genes in these two T-cell groups only:

Seurat::VlnPlot(seu_int, 
            features = c("CD4", "CD8A", "CD8B"),
            idents = c("CD8+ T cells", "CD4+ T cells"))

Differential expression using limma

The Wilcoxon test implemented in FindMarkers does not allow you to test for complex design (eg factorial experiments) or to include batch as a covariate. It doesn’t allow you to run paired-sample T tests for example.

For more complex designs, we can use edgeR or limma which are designed for microarray or bulk RNA seq data and provide a design matrix that includes covariates for example, or sample IDs for paired analyses.

We will load an object containing only pro B cells, both from healthy tissues (PBMMC), and malignant tissues (ETV6-RUNX1).

Warning

Please NOTE that in the original design of this data set, the healthy and malignant tissues were not patient-matched, i.e. the real design was not the one of paired healthy and malignant tissues. However, for demonstration purposes, we will show you how to run a paired analysis, and do as if the PBMMC-1 and ETV6-RUNX1-1 samples both came from the same patient 1, the PBMMC-2 and ETV6-RUNX1-2 samples both came from the same patient 2, etc…

We can load the object and explore its UMAP and meta.data like this:

proB <- readRDS("course_data/proB.rds")

DimPlot(proB, group.by = "orig.ident")

table(proB@meta.data$type)

ETV6-RUNX1      PBMMC 
      2000       1021 
# ETV6-RUNX1      PBMMC 
#      2000       1021

head(proB@meta.data)
                           orig.ident nCount_RNA nFeature_RNA    SingleR_annot
PBMMC-1_AAATGCCAGACTGGGT-1    PBMMC-1       4886         1727 Pro-B_cell_CD34+
PBMMC-1_AAATGCCTCCACTGGG-1    PBMMC-1       8397         2291 Pro-B_cell_CD34+
PBMMC-1_AACACGTTCTTGACGA-1    PBMMC-1       3444         1204 Pro-B_cell_CD34+
PBMMC-1_AACCATGAGAAGGTGA-1    PBMMC-1       8981         2437 Pro-B_cell_CD34+
PBMMC-1_AACCGCGCATGGTCAT-1    PBMMC-1       3719         1368 Pro-B_cell_CD34+
PBMMC-1_AAGCCGCCAGACGTAG-1    PBMMC-1       4573         1464 Pro-B_cell_CD34+
                            type
PBMMC-1_AAATGCCAGACTGGGT-1 PBMMC
PBMMC-1_AAATGCCTCCACTGGG-1 PBMMC
PBMMC-1_AACACGTTCTTGACGA-1 PBMMC
PBMMC-1_AACCATGAGAAGGTGA-1 PBMMC
PBMMC-1_AACCGCGCATGGTCAT-1 PBMMC
PBMMC-1_AAGCCGCCAGACGTAG-1 PBMMC
Note

If you want to know how this pro-B cell subset is generated, have a look at the script here.

Since we will start with differential gene expression, we set the default assay back to “RNA”. Also, we set the default identity to the cell type:

Seurat::DefaultAssay(proB) <- "RNA"
Seurat::Idents(proB) <- proB$orig.ident

Let’s have a look at the UMAP (again), coloured by celltype:

Seurat::DimPlot(proB)

Let’s say we are specifically interested to test for differential gene expression between the tumor and normal samples.

Note

Here we could also test for e.g. healthy versus diseased within a celltype/cluster.

Now we will run differential expression analysis between tumor and healthy cells using the patient ID as a covariate by using limma.

Prepare the pseudobulk count matrix:

#taking the proB data 
Seurat::DefaultAssay(proB) <- "RNA"
Seurat::Idents(proB) <- proB$orig.ident

## add the patient id also for paired DGE
proB$patient.id<-gsub("ETV6-RUNX1", "ETV6_RUNX1", proB$orig.ident)
proB$patient.id<-sapply(strsplit(proB$patient.id, "-"), '[', 2)

## Here we do perform pseudo-bulk:
##first a mandatory column of sample needs to be added to the meta data that is the grouping factor, should be the samples
proB$sample <- factor(proB$orig.ident)

##first an sce object is needed
sce_proB <- as.SingleCellExperiment(proB)

##aggregateAcrossCells here it is only aggregated by sample, one could imagine
##to aggregate by sample and by celltype for instance
summed <- aggregateAcrossCells(sce_proB, 
                               id=colData(sce_proB)[,c("sample")])
##have a look at the counts
counts(summed)[1:3,]
             ETV6-RUNX1-1 ETV6-RUNX1-2 ETV6-RUNX1-3 PBMMC-1 PBMMC-2 PBMMC-3
RP11-34P13.7            0            0            0       2       0       0
FO538757.3              0            0            0       0       0       0
FO538757.2            138          275           74     129      40     112
#have a look at the colData of our new object summed, can you see type and 
#patient.id are there
head(colData(summed))
DataFrame with 6 rows and 10 columns
               orig.ident nCount_RNA nFeature_RNA    SingleR_annot        type
              <character>  <numeric>    <integer>      <character> <character>
ETV6-RUNX1-1 ETV6-RUNX1-1         NA           NA Pro-B_cell_CD34+  ETV6-RUNX1
ETV6-RUNX1-2 ETV6-RUNX1-2         NA           NA Pro-B_cell_CD34+  ETV6-RUNX1
ETV6-RUNX1-3 ETV6-RUNX1-3         NA           NA Pro-B_cell_CD34+  ETV6-RUNX1
PBMMC-1           PBMMC-1         NA           NA Pro-B_cell_CD34+       PBMMC
PBMMC-2           PBMMC-2         NA           NA Pro-B_cell_CD34+       PBMMC
PBMMC-3           PBMMC-3         NA           NA Pro-B_cell_CD34+       PBMMC
              patient.id       sample        ident          ids    ncells
             <character>     <factor>     <factor>     <factor> <integer>
ETV6-RUNX1-1           1 ETV6-RUNX1-1 ETV6-RUNX1-1 ETV6-RUNX1-1       555
ETV6-RUNX1-2           2 ETV6-RUNX1-2 ETV6-RUNX1-2 ETV6-RUNX1-2      1110
ETV6-RUNX1-3           3 ETV6-RUNX1-3 ETV6-RUNX1-3 ETV6-RUNX1-3       335
PBMMC-1                1 PBMMC-1      PBMMC-1      PBMMC-1            461
PBMMC-2                2 PBMMC-2      PBMMC-2      PBMMC-2            153
PBMMC-3                3 PBMMC-3      PBMMC-3      PBMMC-3            407

Generate a DGEList object to use as input for limma and filter the genes to remove lowly expressed genes. How many are left?

#As in the standard limma analysis generate a DGE object

y <- DGEList(counts(summed), samples=colData(summed)$sample)

##filter lowly expressed (recommanded for limma)
keep <- filterByExpr(y, group=summed$type)
y <- y[keep,]

##see how many genes were kept 
summary(keep)
   Mode   FALSE    TRUE 
logical   11086   10017

Generate a design matrix, including patient ID to model for a paired analysis. If you need help to generate a design matrix, check out the very nice edgeR User Guide, sections 3.3 and 3.4. Extract the sample ID from the meta.data, then create the design matrix:

## Create the design matrix and include the technology as a covariate:
design <- model.matrix(~0 + summed$type + summed$patient.id)

# Have a look
design
  summed$typeETV6-RUNX1 summed$typePBMMC summed$patient.id2 summed$patient.id3
1                     1                0                  0                  0
2                     1                0                  1                  0
3                     1                0                  0                  1
4                     0                1                  0                  0
5                     0                1                  1                  0
6                     0                1                  0                  1
attr(,"assign")
[1] 1 1 2 2
attr(,"contrasts")
attr(,"contrasts")$`summed$type`
[1] "contr.treatment"

attr(,"contrasts")$`summed$patient.id`
[1] "contr.treatment"
# change column/rownames names to more simple group names: 
colnames(design) <- make.names(c("ETV6-RUNX1", "PBMMC","patient2","patient3"))
rownames(design)<-colData(summed)$sample

Specify which contrast to analyse:

contrast.mat <- limma::makeContrasts(ETV6.RUNX1 - PBMMC,
                                     levels = design)

Firt, we perform TMM normalization using edgeR, and then limma can perform the transformation with voom, fit the model, compute the contrasts and compute test statistics with eBayes:

dge <- edgeR::calcNormFactors(y)  

#Do limma
vm <- limma::voom(dge, design = design, plot = TRUE)

fit <- limma::lmFit(vm, design = design)
fit.contrasts <- limma::contrasts.fit(fit, contrast.mat)
fit.contrasts <- limma::eBayes(fit.contrasts)

We can use topTable to get the most significantly differentially expressed genes, and save the full DE results to an object. How many genes are significant? Are you suprised by this number?

# Show the top differentially expressed genes:
limma::topTable(fit.contrasts, number = 10, sort.by = "P")
               logFC  AveExpr         t      P.Value    adj.P.Val        B
RPS4Y2      5.346800 6.347826  15.39674 3.361453e-08 0.0001152727 9.477129
SDC2        9.070465 2.708711  15.33434 3.493119e-08 0.0001152727 7.681010
IGLL1      -3.788160 9.287148 -15.19483 3.808426e-08 0.0001152727 9.465422
CTGF        4.368363 6.029640  14.89301 4.603081e-08 0.0001152727 9.141505
AP005530.2  8.770808 2.560369  14.27130 6.879238e-08 0.0001335825 7.294079
GNG11       3.500250 6.457777  13.70183 1.008304e-07 0.0001335825 8.495288
HLA-DQA1    2.982748 7.410939  13.39202 1.249067e-07 0.0001335825 8.287394
PTP4A3      3.734865 5.449894  13.23439 1.395246e-07 0.0001335825 8.126807
CD27        4.115561 5.843582  13.23011 1.399471e-07 0.0001335825 8.147433
ALOX5       4.142010 5.682900  13.16658 1.463814e-07 0.0001335825 8.098199
limma_de <- limma::topTable(fit.contrasts, number = Inf, sort.by = "P")
length(which(limma_de$adj.P.Val<0.05))
[1] 2738

And we can check whether this corresponds to the counts by generating a violin plot, or a gene downregulated in tumor, or a gene upregulated in tumor:

Seurat::VlnPlot(proB, "S100A9", split.by = "type")

Seurat::VlnPlot(proB, "SOCS2", split.by = "type")

We can run a similar analysis with Seurat, but this will not take into account the paired design. Run the code below.

tum_vs_norm <- Seurat::FindMarkers(proB, 
                                   ident.1 = "ETV6-RUNX1", 
                                   ident.2 = "PBMMC", 
                                   group.by = "type")
tum_vs_norm <- subset(tum_vs_norm, tum_vs_norm$p_val_adj<0.05)
Exercise (extra)

How many genes are significant? How does the fold change of these genes compare to the fold change of the top genes found by limma?

dim(tum_vs_norm) 
[1] 1893    5

We find 1893 significant genes. If we merge the FindMarkers and the limma results, keep limma’s most significant genes and plot:

merge_limma_FindMarkers <- merge(tum_vs_norm, limma_de, by="row.names",
                           all.x=T)

par(mar=c(4,4,4,4))
plot(merge_limma_FindMarkers$avg_log2FC,
    merge_limma_FindMarkers$logFC,
    xlab="log2FC Wilcoxon", ylab="log2FC limma",
    pch=15, cex=0.5)
abline(a=0, b=1, col="red")

Keep the object

Keep the tum_vs_norm and limma_de objects because we will use this output later for the enrichment analysis in the next section.